
*********************************************************************************
* (1.0) Load analysis sample 
*********************************************************************************

	use 	"$f_2PKmotherpanel", clear	
			
		keep if rural == 1
		gen svydate		= (88*12)+7
				
	tab educ, gen(edu)
	
	gen ref_year_age = current_age
	drop if ref_year_age<15 | ref_year_age>=50

	cap drop agegroup
	gen agegroup = . 
		replace agegroup = 1 if ref_year_age >= 15 & ref_year_age <= 19
		replace agegroup = 2 if ref_year_age >= 20 & ref_year_age <= 24
		replace agegroup = 3 if ref_year_age >= 25 & ref_year_age <= 29
		replace agegroup = 4 if ref_year_age >= 30 & ref_year_age <= 34
		replace agegroup = 5 if ref_year_age >= 35 
		
	char agegroup[omit] 1

	gen LLF = year >= LLFyear
	cap drop LLFagegroup
	gen LLFagegroup = LLF*agegroup
	tab LLFagegroup, gen(LLFagegroup_)
	
	rename id id_orig
	egen id = group(id)
	
	keep if year >= 1964  & year <1980
	gen years_since_16 = years_since_LLF + 16
	char years_since_16[omit] 15

	merge m:1 prov year using  "$d_data\ProductionData.dta", gen(productionmerge)
	cap drop _merge
	merge m:1 prov year using `"$d_data\IMR.dta"'
	
	cap drop _merge
	merge m:1 prov using `"$d_data\sentdown_mean.dta"'
	drop _merge 
	
	gen 	culturalrev = 0 
	replace culturalrev = mean_sentdown if year >= 1966 & year <=1969
	
	do `"$d_do\LLF_StackData.do"'
	
	gen years_since_LLF = year - LLFyear 

	keep if year >= 1964  & year <1980
	gen years_since_16 = years_since_LLF + 16
	char years_since_16[omit] 15

	forval i = 	1/26 { 
		forval j = 1/5 { 
			gen eventyr_`i'_group_`j' = 0
			replace eventyr_`i'_group_`j' = 1 if years_since_16 == `i' & agegroup == `j'
			} 
		} 
		
	forval i = 	1/26 { 
		drop eventyr_`i'_group_1  
		} 	
	forval j = 	2/5{ 
		drop eventyr_15_group_`j'  
		} 

	tab years_since_16, gen(years_since_fe)
	tab agegroup, gen(agegroup_fe)
		
*********************************************************************************
* (2.0) Regression adjusted life tables measure LLF effect in event time 
*********************************************************************************
			
	drop if year == 1968 
	drop if year == 1967
	
	cap drop included 
	gen included = 0
	
	gen married = 0
	replace married = 1 if year >= year_marriage
	
	 drop if years_since_LLF>7
	 char FE_prov[omit] 108 

	*******************************************************
	* (2.1) Outupt unadjusted TFRs
	*******************************************************
	do  `"$d_do\LLF_compute_unadjustedTFR_eventtime.do"'
	 
	 tempfile main
	 save `main'
		
	*******************************************************
	* (2.2) First births
	*******************************************************
	
	use `main', clear 
	
	global controls1 	"m_year_birth i.educ married U5MR_5yr_ave"
	global controls2 	"grain agric_production gdp pop_total_agric culturalrev"

	xi: logit delivered i.years_since_16 i.agegroup eventyr_* $controls1 $controls2 i.year i.prov if parity==0, cluster(prov) nrtolerance(1e-3)
		replace included = 1 if e(sample) 

	* Compute life table inputs implied by regression output 
	preserve 
	
	keep if e(sample) 
	keep if years_since_LLF <= -1 
	
	* Pre-LLF (event year < -1)
	forval ey = 1/14 {
	
		local eventyear = `ey' - 16
		
		* Set event year indicators for years other than the reference year to 0
		forval i = 1/26 { 
			 cap replace _Iyears_sin_`i' = 0 
			 } 		

		cap replace _Iyears_sin_`ey' = 1 

		* Replace age group * event year indicators for years other than reference year_marriage
		* with zero and predict reference group
		
		forval i = 2/5 { 			
			forval j = 1/14 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 
			forval j = 16/26 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 	
			replace eventyr_`ey'_group_`i' = 1 if agegroup == `i'
			} 
		
		cap drop xb
		predict xb, xb
			forval age = 15/49 { 
				sum xb if current_age == `age' 
				local xb = r(mean)
				if `age' == 15 {
					mat XB_`ey' = `eventyear', `age', `xb' 
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `xb'
					mat XB_`ey' = XB_`ey' \ temp
					}
				}
		
		cap drop qx
		predict qx
			forval age = 15/49 { 
				sum qx if current_age == `age' 
				local qx = r(mean)
				if `age' == 15 {
					mat Qx_`ey' = `eventyear', `age', `qx' 
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `qx'
					mat Qx_`ey' = Qx_`ey' \ temp
					}
				}
		}
	
	* Year prior to implementation (event eyar -1)
	local ey = 15 
		local eventyear = `ey' - 16
		
		* Set event year indicators for years other than the reference year to 0
		forval i = 1/26 { 
			 cap replace _Iyears_sin_`i' = 0 
			 } 		

		* Replace age group * event year indicators for years other than reference year_marriage
		* with zero and predict reference group
		
		forval i = 2/5 { 			
			forval j = 1/14 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 
			forval j = 16/26 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 
		}
		
		cap drop xb
		predict xb, xb
			forval age = 15/49 { 
				sum xb if current_age == `age' 
				local xb = r(mean)
				if `age' == 15 {
					mat XB_`ey' = `eventyear', `age', `xb' 
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `xb'
					mat XB_`ey' = XB_`ey' \ temp
					}
				}
		
		cap drop qx
		predict qx
			forval age = 15/49 { 
				sum qx if current_age == `age' 
				local qx = r(mean)
				if `age' == 15 {
					mat Qx_`ey' = `eventyear', `age', `qx'
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `qx'
					mat Qx_`ey' = Qx_`ey' \ temp
					}
				}	
	
	* LLF period (event year >= 0)
	forval ey = 16/26 {
	
		local eventyear = `ey' - 16
		
		* Set event year indicators for years other than the reference year to 0
		forval i = 1/26 { 
			 cap replace _Iyears_sin_`i' = 0 
			 } 		

		cap replace _Iyears_sin_`ey' = 1 

		* Replace age group * event year indicators for years other than reference year_marriage
		* with zero and predict reference group
		
		forval i = 2/5 { 			
			forval j = 1/14 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 
			forval j = 16/26 { 
				replace eventyr_`j'_group_`i' = 0 	
				} 
				
			replace eventyr_`ey'_group_`i' = 1 if agegroup == `i'
			} 
		
		cap drop xb
		predict xb, xb
			forval age = 15/49 { 
				sum xb if current_age == `age' 
				local xb = r(mean)
				if `age' == 15 {
					mat XB_`ey' = `eventyear', `age', `xb' 
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `xb'
					mat XB_`ey' = XB_`ey' \ temp
					}
				}
		
		cap drop qx
		predict qx
			forval age = 15/49 { 
				sum qx if current_age == `age' 
				local qx = r(mean)
				if `age' == 15 {
					mat Qx_`ey' = `eventyear', `age', `qx' 
					} 
				else if `age' > 15 { 
					mat temp = `eventyear', `age', `qx'
					mat Qx_`ey' = Qx_`ey' \ temp
					}
				}

		}
	restore	
	
	forval ey = 11/24 { 
		forval row = 1/35 { 
			di "EY `ey' Row `row'"
			if Qx_`ey'[`row', 3] == . { 
				matrix Qx_`ey'[`row', 3] = 0 
				}
			}
		}
		
	forval i = 11/24 { 	
		local eventyear = `i' - 16 
		mat Lx_`i' = `eventyear', 15, 1

		scalar dx = Qx_`i'[1, 3]*Lx_`i'[1, 3] 
		mat Dx_`i' = `eventyear', 15, dx 
	
		forval j = 2/35 { 
			local age = `j' + 14
			local n_1 = `j' - 1
	
			scalar lx = Lx_`i'[`n_1', 3] - Dx_`i'[`n_1', 3] 
			mat temp = `eventyear', `age', lx
			mat Lx_`i' = Lx_`i' \ temp
	
			scalar dx = Qx_`i'[`j', 3]*Lx_`i'[`j', 3] 
			mat temp = `eventyear', `age' , dx 
			mat Dx_`i' = Dx_`i' \ temp	
			}	
		}
 	
	*******************************************************
	* (2.2) Higher Order births
	*******************************************************

	forval parity = 1/6 { 

		xi: logit delivered i.years_since_16 i.agegroup eventyr_* $controls1 $controls2 i.year i.prov if parity==`parity', cluster(prov)  iter(30) nrtolerance(1e-3)
		replace included = 1 if e(sample) 
		*outreg2 using "$d_tab\Lifetable_regoutput.xls", excel  stats(coef  ) bdec(4) noaster append	

		preserve 
		
		keep if e(sample)
		keep if years_since_LLF <= -1 
			
		forval ey = 1/14 {
		
			local eventyear = `ey' - 16
				
			* Set event year indicators for years other than the reference year to 0
			forval i = 1/26 { 
				 cap replace _Iyears_sin_`i' = 0 
				 } 		

			cap replace _Iyears_sin_`ey' = 1 


			* Replace age group * event year indicators for years other than reference year_marriage
			* with zero and predict reference group
			forval i = 2/5 { 			
				forval j = 1/14 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
				forval j = 16/26 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
					
				replace eventyr_`ey'_group_`i' = 1 if agegroup == `i'
				} 
			
			cap drop xb
			predict xb, xb
				forval age = 15/49 { 
					sum xb if current_age == `age' 
					local xb = r(mean)
					if `age' == 15 {
						mat tXB_`ey' = `xb' 
						} 
					else if `age' > 15 { 
						mat temp = `xb'
						mat tXB_`ey' = tXB_`ey' \ temp
						}
					}
			
			cap drop qx
			predict qx
				forval age = 15/49 { 
					sum qx if current_age == `age' 
					local qx = r(mean)
					if `age' == 15 {
						mat tQx_`ey' = `qx' 
						} 
					else if `age' > 15 { 
						mat temp = `qx'
						mat tQx_`ey' = tQx_`ey' \ temp
						}
					}
		
			mat Qx_`ey' = Qx_`ey', tQx_`ey'
			mat XB_`ey' = XB_`ey', tXB_`ey'

			}	
		
		local ey = 15 
		
			local eventyear = `ey' - 16

			* Set event year indicators for years other than the reference year to 0
			forval i = 1/26 { 
				 cap replace _Iyears_sin_`i' = 0 
				 } 		

			* Replace age group * event year indicators for years other than reference year_marriage
			* with zero and predict reference group
			forval i = 2/5 { 			
				forval j = 1/14 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
				forval j = 16/26 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
			}
			
			cap drop xb
			predict xb, xb
				forval age = 15/49 { 
					sum xb if current_age == `age' 
					local xb = r(mean)
					if `age' == 15 {
						mat tXB_`ey' = `xb' 
						} 
					else if `age' > 15 { 
						mat temp = `xb'
						mat tXB_`ey' = tXB_`ey' \ temp
						}
					}
			
			cap drop qx
			predict qx
				forval age = 15/49 { 
					sum qx if current_age == `age' 
					local qx = r(mean)
					if `age' == 15 {
						mat tQx_`ey' = `qx' 
						} 
					else if `age' > 15 { 
						mat temp = `qx'
						mat tQx_`ey' = tQx_`ey' \ temp
						}
					}
		
			mat Qx_`ey' = Qx_`ey', tQx_`ey'
			mat XB_`ey' = XB_`ey', tXB_`ey'		
		
		
		forval ey = 16/26 {
		
			local eventyear = `ey' - 16
			
			* Set event year indicators for years other than the reference year to 0
			forval i = 1/26 { 
				 cap replace _Iyears_sin_`i' = 0 
				 } 		

			cap replace _Iyears_sin_`ey' = 1 

			* Replace age group * event year indicators for years other than reference year_marriage
			* with zero and predict reference group
			forval i = 2/5 { 			
				forval j = 1/14 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
				forval j = 16/26 { 
					replace eventyr_`j'_group_`i' = 0 	
					} 
					
				replace eventyr_`ey'_group_`i' = 1 if agegroup == `i'
				} 
			
			cap drop xb
			predict xb, xb
				forval age = 15/49 { 
					sum xb if current_age == `age' 
					local xb = r(mean)
					if `age' == 15 {
						mat tXB_`ey' = `xb' 
						} 
					else if `age' > 15 { 
						mat temp = `xb'
						mat tXB_`ey' = tXB_`ey' \ temp
						}
					}
			
			cap drop qx
			predict qx
				forval age = 15/49 { 
					sum qx if current_age == `age' 
					local qx = r(mean)
					if `age' == 15 {
						mat tQx_`ey' = `qx' 
						} 
					else if `age' > 15 { 
						mat temp = `qx'
						mat tQx_`ey' = tQx_`ey' \ temp
						}
					}
		
			mat Qx_`ey' = Qx_`ey', tQx_`ey'
			mat XB_`ey' = XB_`ey', tXB_`ey'

			}
		restore	
		}
		
		forval parity = 1/6 { 
		forval ey = 11/24 { 
			forval row = 1/35 { 
			local col_p = `parity' + 3
				if Qx_`ey'[`row', `col_p'] == . { 
					matrix Qx_`ey'[`row', `col_p'] = 0 
					}
				}
			}
		}
		
		forval parity = 1/6 { 
		forval i = 11/24 { 
		
			local col_p = `parity' + 3
			local col_p_1 = `col_p' - 1
		
			local r = rowsof(Lx_`i')
			mat temp = J(`r', 1, .)
		
			mat Lx_`i' = Lx_`i', temp 
			mat Dx_`i' = Dx_`i', temp 
		
			mat Lx_`i'[1, `col_p'] =  0
			mat Dx_`i'[1, `col_p'] = Qx_`i'[1, `col_p']*(Lx_`i'[1, `col_p']+(Dx_`i'[1, `col_p_1']/2)) 
		
			forval n = 2/35 { 
				local n_1 = `n' - 1
		
				mat Lx_`i'[`n', `col_p'] = Lx_`i'[`n_1', `col_p'] + Dx_`i'[`n_1', `col_p_1'] - Dx_`i'[`n_1', `col_p'] 
				mat Dx_`i'[`n', `col_p'] = Qx_`i'[`n', `col_p']*(Lx_`i'[`n', `col_p']+(Dx_`i'[`n', `col_p_1']/2))	
				}	
			}
		}

*******************************************************
* (3.0) Output key life table components
*******************************************************

	forval i = 11/24 { 
		mat colnames Qx_`i' = "Eventyear" "Age" "Qx_1" "Qx_2" "Qx_3" "Qx_4" "Qx_5" "Qx_6" "Qx_7"
		mat colnames Lx_`i' = "Eventyear" "Age" "Lx_1" "Lx_2" "Lx_3" "Lx_4" "Lx_5" "Lx_6" "Lx_7"
		mat colnames Dx_`i' = "Eventyear" "Age" "Dx_1" "Dx_2" "Dx_3" "Dx_4" "Dx_5" "Dx_6" "Dx_7"
		}
				
	cap mkdir `"$d_data\Fertility Life Tables"' 
	
	preserve 		
	forval i = 11/24 { 
		clear 
		svmat Qx_`i', names(col) 
			save `"$d_data\Fertility Life Tables\Qx_`i'.dta"', replace
		clear 
		svmat Lx_`i', names(col) 
			save `"$d_data\Fertility Life Tables\Lx_`i'.dta"', replace
		clear 
		svmat Dx_`i', names(col) 
			gen ASFR = Dx_1 + Dx_2 + Dx_3 + Dx_4 + Dx_5 + Dx_6 + Dx_7
			save `"$d_data\Fertility Life Tables\Dx_`i'.dta"', replace

		clear 
		use `"$d_data\Fertility Life Tables\Dx_`i'.dta"'
				
			count 
			local max = r(N) 
			
			scalar TFR_`i' = ASFR[1]
			forval parity = 1/7 { 
				scalar PPR`parity'_`i' = Dx_`parity'[1]
				}
			forval age = 2/35 { 
				scalar TFR_`i' = TFR_`i' + ASFR[`age']
				forval parity = 1/7 { 
					scalar PPR`parity'_`i' = PPR`parity'_`i' + Dx_`parity'[`age']
					} 
				}
					
			if `i' == 11 {
				mat TFR = TFR_`i'
				
				mat PPR = 1, PPR1_`i'
				forval parity = 2/7 { 
					mat temp = `parity', PPR`parity'_`i' 
					mat PPR = PPR \ temp
					}
				} 
				
			else if `i' > 11 { 
				mat TFR = TFR, TFR_`i' 
				
				mat temp1 = 1, PPR1_`i'
				forval parity = 2/7 { 
					mat temp2 = `parity', PPR`parity'_`i' 
					mat temp1 = temp1 \ temp2
					}
				mat PPR = PPR, temp1
				} 
			
		* Output age specific fertility rate (ASFR)
		keep Age ASFR
		rename ASFR ASFR_`i' 
		if `i' == 11 {
			save `"$d_data\Fertility Life Tables\ASFRs.dta"', replace
			}
		else if `i' > 11 { 
			merge 1:1 Age using `"$d_data\Fertility Life Tables\ASFRs.dta"', nogen
			save `"$d_data\Fertility Life Tables\ASFRs.dta"', replace
			}
		}
	
*******************************************************
* (4.0) Output TFRs in event time
*******************************************************
	
	clear 
	svmat TFR 
		rename TFR14 TFR24
		rename TFR13 TFR23
		rename TFR12 TFR22
		rename TFR11 TFR21
		rename TFR10 TFR20
		rename TFR9 TFR19
		rename TFR8 TFR18
		rename TFR7 TFR17
		rename TFR6 TFR16
		rename TFR5 TFR15
		rename TFR4 TFR14
		rename TFR3 TFR13
		rename TFR2 TFR12
		rename TFR1 TFR11 
		
		forval i = 11/24 { 
			local j = `i' -16
			lab var TFR`i' "TFR Event Year `j'"
			} 
			
		save `"$d_data\Fertility Life Tables\TFRs.dta"', replace
	
*******************************************************
* (5.0) Make Tables and Figure
*******************************************************
	*******************************************************
	* (5.1) Figure 3 
	*******************************************************
	 
	use `"$d_data\Fertility Life Tables\TFRs.dta"', clear
	gen i = 1
	reshape long TFR, i(i) j(eventtime)
	gen temp = _n 
	gen eventyear = temp -6
	
	merge 1:1 eventyear using `"$d_data\TFR_undadjusted_eventtime.dta"'
	rename tfrunadjusted TFRobserved
	drop i 
	drop temp
	drop _merge
	
	* TFR in year prior to implementation 
	gen temp = TFR if eventyear == -1
	egen TFRrefyear = mean(temp)
	drop temp
	
	* LLF-related part of the change in TFR (negative)
	gen TFRchange = TFR - TFRrefyear
	
	gen temp = TFRobserved if eventyear == -1
	egen TFRref = mean(temp)
	drop temp
	
	* The implied TFR if there had not been any LLF is the undadjusted TFR - the LLF related change 
	gen impliedTFRnoLLF = TFRobserved - TFRchange
	
	
	replace TFR = round(TFR, .001)
	sum TFR if eventyear == 7
	local post = r(mean)
	local post = round(`post', .001)
	sum TFR if eventyear == 0 
	local pre = r(mean)
	local pre = round(`pre', .001)
	
	sum TFRchange if eventyear==7
	global change = r(mean)
	global change = round($change , .01)
	local pctchange = -($change/3.09)*100
	local pctchange = round(`pctchange', 1)
	
	twoway (connected TFRobserved eventyear, msymbol(dot) msize(zero) lcolor(black) lwidth(medthick)) ///
	(connected impliedTFRnoLLF eventyear if eventyear>=0, msymbol(dot) msize(zero) lwidth(medthick) lpattern(dash) ) if eventyear<8, ///
	yscale(range(3 6)) ylab(3(.25)6, format(%9.2f)) xlab(-5(1)7) xline(0) ///
	xtitle("Event Year") ytitle("Total Fertility Rate") legend(pos(6) order (1 "Observed TFR" 2 "Implied TFR without LLF") rows(1)) ///
	note("Change in TFR associated with LLF: $change births per woman, or `pctchange'% of overall decline")
	graph export `"$d_fig\Fig_3a_TFRchange.jpg"', replace
		
	restore	
	
	* Convert to caldenar year 
	* To do so, we first map LLF-driven TFR changes into calendar year relevant for each province. 
	* Then we population weight those province years and sum to get a calendar year interpretation 
	* of the LLF-driven TFR decline. 

	preserve 
	keep prov LLFyear year years_since_LLF pop_total
	duplicates drop 
	drop if year < 1968 | year > 1978
	tab year 
	
	drop if prov==15 | prov== 46 | prov==33 | prov==11 
	tempfile pop 
	save `pop', replace 
	
	use `"$d_data\Fertility Life Tables\TFRs.dta"', clear 
	xpose, clear 
	gen years_since_LLF=_n-6
	rename v1 TFR
	
	merge 1:m years_since_LLF using `pop', nogen
	bysort year: egen tot_pop = total(pop_total)
	gen prov_yr_wt = pop_total/tot_pop
	gen temp=TFR*prov_yr_wt
	bysort year: egen TFR_LLFcy = total(temp)
	
	keep year TFR_LLFcy
	duplicates drop 

	tempfile llf_tfr
	save `llf_tfr', replace 
		import delimited "$d_data\UNPD_WFD_2017_TFR.csv", clear 
		drop if year < 1967 | year>1978
		merge 1:1 year using `llf_tfr'
		rename tfr True_TFR
		sum TFR_LLFcy if year==1970
			local base = r(mean)
		gen cum_change = TFR_LLFcy-`base'
		gen TFR_woLLF = True_TFR if year<=1970
		replace TFR_woLLF = True_TFR-cum_change if year>1970

	sum cum_change if year==1978 
		local LLFchange = r(mean)
	sum True_TFR if year==1970
		local tfr1970 = r(mean)
	sum True_TFR if year==1978
		local tfr1978 = r(mean)
		local truechange=`trf1970'-`tfr1978'
		local pctchange = (`LLFchange'/`truechange')*100
		local pctchange = round(`pctchange', .1)
		local LLFchange = round(`LLFchange', .001)
	twoway (line True_TFR year if year <1979) (line TFR_woLLF year if year <1979 & year>=1970, lpattern(dash)), ///
	legend(order(1 "Observed TFR" 2 "Implied TFR without LLF") pos(6) rows(1)) title("") ///
	xtitle("Calendar Year") ytitle("Total Fertility Rate") xlab(1968(1)1978) xscale(range(1968 1978)) ///
	xline(1970, lcolor(black))  ylab(3(.25)6, format(%9.2f)) note("Change in TFR associated with LLF: `LLFchange' births or `pctchange'% of overall decline")
		
	graph export "$d_fig\Fig_3b_TFRchange.jpg", replace
	

	*******************************************************
	* (5.3) Appendix Figure 8 
	*******************************************************
	use `"$d_data\Fertility Life Tables\Dx_16.dta"', clear 	
		drop ASFR
		drop Eventyear
		forval i = 1/7 {
			rename Dx_`i' Dx_`i'_0
		}
		
	merge 1:1 Age using `"$d_data\Fertility Life Tables\Dx_23.dta"'
		drop ASFR
		drop Eventyear
		drop _merge
		forval i = 1/7 {
			rename Dx_`i' Dx_`i'_7
		}		
	
	forval i = 1/7 {
		gen dx_change_`i' =  Dx_`i'_7 - Dx_`i'_0 
	}

	keep Age dx_change_*
	reshape long dx_change_, i(Age) j(parity)	
	keep if parity==1
	gen dx_change_smoothed = . 
	replace dx_change_smoothed = (dx_change_+dx_change_[_n+1])/2 if Age==15
	replace dx_change_smoothed = (dx_change_[_n-1]+dx_change_+dx_change_[_n+1])/3 if Age>15 & Age<=45
	
	twoway (connected dx_change_smoothed Age if parity==1, mcolor(midblue) msize(tiny) ///
	msymbol(point) lcolor(midblue)) if Age>=15 & Age<=40, ///
	xlab(15(1)40, angle(45)) yline(0) ylab(-.015(.005).02, format(%05.3f)) ///
	ytitle("") xtitle("") title("Change in Age Specific First Parity Fertility Rate" "Comparing Year of LLF Implementation to Seven Years Post-LLF", size(msmall))
	
	graph export "$d_fig\Fig_A8_FirstParityASFR.jpg", replace	

	
	*******************************************************
	* (5.4) Appendix Table 8
	*******************************************************

	putexcel set  `"$d_tab\Table_A8.xlsx"', replace 
	
		putexcel B1:F1 = "Age Specific Fertility Rates", merge border(bottom) hcenter
		putexcel H1:I1 = "Total Fertility Rate", merge border(bottom) hcenter
		putexcel B2 = "15-19", hcenter
		putexcel C2 = "20-24", hcenter
		putexcel D2 = "25-29", hcenter
		putexcel E2 = "30-34", hcenter
		putexcel F2 = "35+", hcenter
		putexcel H2 = "TFR", hcenter
		putexcel I2 = "TFR'", hcenter
		putexcel A2:F2, border(bottom)
		putexcel H2:I2, border(bottom)
		putexcel A3 = "5 Years Prior To LLF"
		putexcel A4 = "4 Years Prior To LLF"
		putexcel A5 = "3 Years Prior To LLF"
		putexcel A6 = "2 Years Prior To LLF"
		putexcel A7 = "1 Year Prior To LLF"
		putexcel A8 = "Year of LLF"
		putexcel A9 = "1 Year After LLF"
		putexcel A10 = "2 Years After LLF"
		putexcel A11 = "3 Years After LLF"
		putexcel A12 = "4 Years After LLF"
		putexcel A13 = "5 Years After LLF"
		putexcel A14 = "6 Years After LLFF"
		putexcel A15 = "7 Years After LLF"
		putexcel A15:I15, border(bottom)

	* ASFRs Stored here 
	use `"$d_data\Fertility Life Tables\ASFRs.dta"', clear 	
	gen agegroup = 1 if Age<=19
		replace agegroup = 2 if Age>19 & Age<=24
		replace agegroup = 3 if Age>24 & Age<=29
		replace agegroup = 4 if Age>29 & Age<=34
		replace agegroup = 5 if Age>34
	drop Age
	
	collapse (sum) ASFR_11 ASFR_12 ASFR_13 ASFR_14 ASFR_15 ASFR_16 ///
	ASFR_17 ASFR_18 ASFR_19 ASFR_20 ASFR_21 ASFR_22 ASFR_23 ASFR_24, by(agegroup)
	
	xpose, clear  
	drop in 1 
	drop in 14
	forval i = 1/5 {
		replace v`i' = round(v`i', .01)
	}
	mkmat v1 v2 v3 v4 v5, mat(ASFR) 
	putexcel B3 = mat(ASFR)
	
	use `"$d_data\Fertility Life Tables\TFRs.dta"', clear
	xpose, clear  
	drop in 14	
	replace v1 = round(v1, .01)
	mkmat v1, mat(TFR)
	putexcel H3 = mat(TFR) 
	
	* Make the parity specific TFR (PSFR) so we can compute Bongaarts Feeney decomposition
	use `"$d_data\Fertility Life Tables\Dx_11.dta"', clear
	
	forval i = 1/7 { 
		cap drop P`i'FR
		egen P`i'FR = total(Dx_`i')
		} 
	keep Eventyear P1FR P2FR P3FR P4FR P5FR P6FR P7FR
		duplicates drop 
		save `"$d_data\Fertility Life Tables\PSFRs.dta"', replace 	
	
	forval yr = 12/24 { 
		use `"$d_data\Fertility Life Tables\Dx_`yr'.dta"'
	
		forval i = 1/7 { 
			cap drop P`i'FR
			egen P`i'FR = total(Dx_`i')
			} 
		keep Eventyear P1FR P2FR P3FR P4FR P5FR P6FR P7FR
		duplicates drop 
		
		append using `"$d_data\Fertility Life Tables\PSFRs.dta"'
		save `"$d_data\Fertility Life Tables\PSFRs.dta"', replace 	

		}
				
		use `"$d_data\Fertility Life Tables\PSFRs.dta"', clear 
		sort Eventyear
		
	* Divide by (1 - ave change in age of first birth) 
		gen P1FR_BF = P1FR/(1-0.0366)
		gen TFR_BF = P1FR_BF + P2FR + P3FR + P4FR + P5FR + P6FR + P7FR
		keep if Eventyear>=0 & Eventyear<=7
		keep TFR_BF
		replace TFR_BF = round(TFR_BF, .01)
		mkmat TFR_BF, mat(TFR_BF_post)
		mat TFR_BF=TFR[1..5, 1]
		mat TFR_BF=TFR_BF \ TFR_BF_post
		putexcel I3 = mat(TFR_BF) 
		
	putexcel save 

	*******************************************************
	* (5.5) Appendix Table 7
	*******************************************************
	
	putexcel set `"$d_tab\Table_A7.xlsx"', replace 
		putexcel B1:H1 = "Parity Specific Fertility" , merge border(bottom) hcenter
		putexcel J1:P1 = "Tempo-Adjusted Parity Specific Fertility" , merge border(bottom) hcenter
		putexcel R1:S1 = "Total Fertility Rate" , merge border(bottom) hcenter
		putexcel A2 = "Event"
		putexcel A3 = "Year"
		putexcel B3 = "Parity 1"
		putexcel C3 = "Parity 2"
		putexcel D3 = "Parity 3"
		putexcel E3 = "Parity 4"
		putexcel F3 = "Parity 5"
		putexcel G3 = "Parity 6"
		putexcel H3 = "Parity 7"
		putexcel J3 = "Parity 1"
		putexcel K3 = "Parity 2"
		putexcel L3 = "Parity 3"
		putexcel M3 = "Parity 4"
		putexcel N3 = "Parity 5"
		putexcel O3 = "Parity 6"
		putexcel P3 = "Parity 7"	
		putexcel R3 = "TFR"
		putexcel S3 = "TFR'"
		putexcel A3:S3, border(bottom)

		use `"$d_data\Fertility Life Tables\TFRs.dta"', clear
		xpose, clear  
		drop in 14	
		replace v1 = round(v1, .001)
		mkmat v1, mat(TFR)
		putexcel R4 = mat(TFR)
		
		use `"$d_data\Fertility Life Tables\PSFRs.dta"', clear 
		sort Eventyear
		drop if Eventyear>7
		
		* Parity-specific fertility rates 
		forval i = 1/7 {
			replace P`i'FR = round(P`i'FR, .001)
		}
		
		mkmat Eventyear P1FR P2FR P3FR P4FR P5FR P6FR P7FR, mat(psfr)
			putexcel A4 = mat(psfr)
		
		* Tempo adjust LLF era parity-specific fertility rates (only age at parity 1 birth changes)
		forval i = 1/7 {
			gen P`i'FR_BF = P`i'FR
		}
		replace P1FR_BF = P1FR/(1-0.0366) if Eventyear>=0
		replace P1FR_BF = round(P1FR_BF, .001)

		mkmat P1FR_BF P2FR_BF P3FR_BF P4FR_BF P5FR_BF P6FR_BF P7FR_BF, mat(psfr_bf)
			putexcel J4 = mat(psfr_bf)
			
		gen TFR_BF = P1FR_BF + P2FR + P3FR + P4FR + P5FR + P6FR + P7FR
		keep if Eventyear>=0 
		replace TFR_BF = round(TFR_BF, .001)
		mkmat TFR_BF, mat(TFR_BF_post)
		mat TFR_BF=TFR[1..5, 1]
		mat TFR_BF=TFR_BF \ TFR_BF_post
		putexcel S4 = mat(TFR_BF)
		putexcel A16:S16, border(bottom)
		
	putexcel save
	
	*******************************************************
	* (5.6) Appendix Table 6, all panels
	*******************************************************
	forval i = 1/13 { 
	local ey =`i'-6
	local j = `i' + 10
	putexcel set  `"$d_tab\Table_A6_`i'.xlsx"', replace 
		putexcel A1:Z1 = "Event Year `ey'", merge hcenter bold
		putexcel B2:H2 = "Hazard Function qx", merge border(bottom) hcenter
		putexcel J2:P2 = "Survival Function lx", merge border(bottom) hcenter
		putexcel R2:X2 = "Probability of a Birth dx", merge border(bottom) hcenter
		putexcel Z2 = "Fertility Rate"
		putexcel A3 = "Age"
		putexcel B3 = "Parity 1"
		putexcel C3 = "Parity 2"
		putexcel D3 = "Parity 3"
		putexcel E3 = "Parity 4"
		putexcel F3 = "Parity 5"
		putexcel G3 = "Parity 6"
		putexcel H3 = "Parity 7"
		putexcel J3 = "Parity 1"
		putexcel K3 = "Parity 2"
		putexcel L3 = "Parity 3"
		putexcel M3 = "Parity 4"
		putexcel N3 = "Parity 5"
		putexcel O3 = "Parity 6"
		putexcel P3 = "Parity 7"	
		putexcel R3 = "Parity 1"
		putexcel S3 = "Parity 2"
		putexcel T3 = "Parity 3"
		putexcel U3 = "Parity 4"
		putexcel V3 = "Parity 5"
		putexcel W3 = "Parity 6"
		putexcel X3 = "Parity 7"
		putexcel A3:Z3, border(bottom)
		putexcel A38:Z38, border(bottom)
		
	use `"$d_data\Fertility Life Tables\Qx_`j'.dta"', clear 	
	drop Eventyear
	forval i = 1/7 {
		replace Qx_`i' = round(Qx_`i', .001)
	}
	mkmat Age Qx_1 Qx_2 Qx_3 Qx_4 Qx_5 Qx_6 Qx_7, mat(qx)
	putexcel A4 = mat(qx)

	use `"$d_data\Fertility Life Tables\Lx_`j'.dta"', clear 	
	drop Eventyear Age
	forval i = 1/7 {
		replace Lx_`i' = round(Lx_`i', .001)
	}
	mkmat  Lx_1 Lx_2 Lx_3 Lx_4 Lx_5 Lx_6 Lx_7, mat(lx)
	putexcel J4 = mat(lx)
	
	use `"$d_data\Fertility Life Tables\Dx_`j'.dta"', clear 	
	drop Eventyear Age
	forval i = 1/7 {
		replace Dx_`i' = round(Dx_`i', .001)
	}
	mkmat  Dx_1 Dx_2 Dx_3 Dx_4 Dx_5 Dx_6 Dx_7, mat(dx)
	putexcel R4 = mat(dx)
	
	use `"$d_data\Fertility Life Tables\ASFRs.dta"', clear
	keep ASFR_`j'
	replace ASFR_`j'=round(ASFR_`j', .001)
	mkmat ASFR_`j', mat(asfr)
	putexcel Z4 = mat(asfr)
	
	use `"$d_data\Fertility Life Tables\PSFRs.dta"', clear
	keep if Eventyear==`ey'
	drop Eventyear
	forval i = 1/7 { 
		replace P`i'FR = round(P`i'FR, .001)
	} 
	mkmat P*, mat(psfr)
	putexcel R39 = mat(psfr)
	
	gen tfr = P1FR+P2FR+P3FR+P4FR+P5FR+P6FR+P7FR
	replace tfr=round(tfr,.001)
	sum tfr
	putexcel Z39 = `r(mean)'
	}
	
	putexcel save  	
	
	*******************************************************
	* (5.2) Appendix Figure 16 
	*******************************************************
	use `"$d_data\Fertility Life Tables\Dx_16.dta"', clear 	
		drop ASFR
		drop Eventyear
		forval i = 1/7 {
			rename Dx_`i' Dx_`i'_0
		}
		
	merge 1:1 Age using `"$d_data\Fertility Life Tables\Dx_23.dta"'
		drop ASFR
		drop Eventyear
		drop _merge
		forval i = 1/7 {
			rename Dx_`i' Dx_`i'_7
		}		
	
	forval i = 1/7 {
		gen dx_change_`i' =  Dx_`i'_7 - Dx_`i'_0 
	}

	keep Age dx_change_*
	reshape long dx_change_, i(Age) j(parity)

	* This figure is finicky - Stata doesn't seem to like the RBG color specification. 
	* If it doesn't output correctly, I've found deleting and replacing RBG codes to work 
	* - or you can set your own color scheme 
	
	cap twoway contour dx_change_ Age parity, ccuts(-.05(.01).015) ///
	lcolor(blue) graphregion(fcolor(white) lcolor(black)) ylabel(15(5)50, nogrid angle(horizontal))  ///
	ccolors("0 0 0" "51 2 4" "99 3 6" "113 3 7" "131 3 8" "158 4 10" "197 5 12" "212 68 73" "226 130 134" "255 203 203" "255 255 255" "231 231 253" "187 187 225") ///
	xlab(1(1)7) scale(.75) aspectratio(.75) xtitle("Parity") ccuts(.015, .005, -.005, -.015, -.02, -.025, -.03, -0.035, -.04, -.045, -.05,  -.055) ///
	ytitle("Maternal Age") ztitle("Change in Age and Parity Specific Fertility Rate" "(Births per Woman)") xscale(range(1(1)7)) interp(shepard) crule(linear) scolor(black) ecolor(white)

	gr_edit .clegend.zaxis.major.label_format = `"%9.3f"'
	
	graph export "$d_fig\Fig_A16_ASFRContourMap.jpg", replace	
	
